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Salt transport in bulk electrolytes is limited by diffusion and convection, but in microstructures 
with charged surfaces (e.g. microfluidic devices, porous media, soils, or biological tissues) surface 
conduction and electro-osmotic flow also contribute to ionic fluxes. For small applied voltages, these 
effects lead to well known linear electrokinetic phenomena. In this paper, we predict some surprising 
nonlinear dynamics that can result from the competition between bulk and interfacial transport at 
higher voltages. When counter-ions are selectively removed by a membrane or electrode, a "desali- 
nation shock" can propagate through the microstructure, leaving in its wake an ultrapure solution, 
nearly devoid of co-ions and colloidal impurities. We elucidate the basic physics of desalination 
shocks and develop a mathematical theory of their existence, structure, and stability, allowing for 
slow variations in surface charge or channel geometry. Via asymptotic approximations and similar- 
ity solutions, we show that desalination shocks accelerate and sharpen in narrowing channels, while 
they decelerate and weaken, and sometimes disappear, in widening channels. These phenomena 
may find applications in separations (desalination, decontamination, biological assays) and energy 
storage (batteries, supercapacitors) involving electrolytes in microstructures. 
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I. INTRODUCTION 

All electrochemical processes lead to ionic concentra- 
tion gradients in electrolytes [l|, [2| . In water desalina- 
tion, the removal of ions is the desired outcome, but in 
most other situations, such as energy storage by batter- 
ies or energy conversion by fuel cells, salt depletion is 
undesirable because it increases the solution resistance 
and slows electrochemical reactions, thereby increasing 
the over-potential required to maintain a desired current. 
Salinity variations also commonly arise in biological sys- 
tems due to the action of membranes or external stimuli, 
and their dynamics can significantly affect living cells 
and tissues. In all of these situations it is important to 
understand the dynamics of ions in complex geometries. 

It is generally assumed that salt transport in bulk 
electrolytes occurs only by diffusion and convection. 
This hypothesis underlies important industrial pro- 
cesses, such as electrodialysis [3, 4], electrodeposition [5], 
and experimental techniques, such as impedance spec- 
troscopy [6], cyclic voltammetry [7j. In a concen- 
trated electrolyte, ionic diffusion is nonlinear (with a 
concentration-dependent diffusivity Q), but the famil- 
iar square-root of time scaling of linear diffusion usually 
remains |8j. This conclusion also holds for macroscopic 
transport in porous media, as long as linear diffusion oc- 
curs within the pores [9]. 

Recent experiments have shown that more compli- 



cated, nonlinear dynamics are possible if strong salt 
depletion ("concentration polarization") occurs in mi- 
crostructures. A growing body of work has fo- 
cused on Dukhin's second-kind electro-osmotic flows 
fllj and the Rubinstein- Zaltzman instability [HI, 
near electrodialysis membranes 0, [HJ and microchan- 
nel/nano channel junctions [HI, [HI and in packed beds 
of particles [13, Ho]- ^ n a ^ °f these cases, the transport 
of ions across a selective surface depletes the salt concen- 
tration and causes nonlinear electrokinetic phenomena in 
electric double layers (EDLs) sustaining normal current. 

In contrast, our focus here is on the effect of tangen- 
tial current in the EDL [l9l-[2lj. also known as "surface 
conduction" , which has a long history, prior to microflu- 
idics [22l-[2(|. In linear electrokinetics, the importance 
of surface conduction is controlled by the Dukhin num- 
ber M Szi , 
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where is the conductivity of the neutral bulk solution, 
k' s is the additional "surface conductivity" due to excess 
ions in the EDLs [H [H, 0, HI , and h is a geometrical 
length scale, such as the channel width or particle size. 
The competition of surface and bulk conduction in a mi- 
crochannel is now well understood for linear response to 
a small voltage or current [H, [28|, [29|, but recently a 
surprising nonlinear phenomenon was discovered. 

Mani, Zangle and Santiago showed that, under certain 
conditions, surface conduction can produce a localized 
salt concentration gradient propagating through a mi- 
crochannel, away from a nanochannel junction (3(J l3lj |. 
By deriving a one-dimensional equation for thin EDLs 
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(the "Simple Model" ) and applying the method of char- 
acteristics, they explained this phenomenon mathemat- 
ically as shock propagation in the concentration profile, 
analogous to pressure shocks in gases [30|. The theory 
was able to predict, for the first time, the propagation 
of enrichment and depletion shocks in etched glass mi- 
crochannels on either side of a nanochannel [3l|. It is 
possible that this phenomenon plays a role in earlier ob- 
servations of sharp concentration gradients in more com- 
plicated microchannel/nanochannel geometries [HI, 32- 

In this paper, we focus on the new surface-conduction 
dominated regime and develop a general theory of "de- 
salination shocks" in complex microstructures. We be- 
gin by describing the basic physics of desalination shock 
propagation in microchannels or porous media. We then 
develop general macroscopic transport equations for ions 
in charged microstructures, which lead to a nonlinear 
wave equation at constant current. After making the 
equations dimensionless and identifying the key govern- 
ing parameters, we study desalination shock propagation 
in two types of heterogeneous microstructures. First, we 
analyze slowly varying surface charge and/or channel ge- 
ometry using perturbation methods, and then we derive 
intermediate-asymptotic similarity solutions for power- 
law variations in the channel area. The latter clarify the 
transition from diffusive scaling (x ~ y/t) without shocks 
in a wedge to constant-velocity shock propagation in a 
straight channel (x ~ t). Finally, we show that thin de- 
salination shocks are nonlinearly stable in the absence of 
fluid flow by reducing the dynamics to a Laplacian disso- 
lution model. We conclude by discussing possible appli- 
cations of our results to microfluidic separations, water 
desalination, soil decontamination, and energy storage by 
porous electrodes. 
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FIG. 1: Basic physics of desalination shock propagation. (A) 
Sketch of ion fluxes in a microchannel or pore with negatively 
charged walls, as current flows from left to right through a de- 
crease in salt concentration (caused by an electrode or mem- 
brane, not shown). In order to avoid low- conductivity region 
in the center of the channel, the current flows into the elec- 
tric double layers, where it is carried by positive counter- ions 
that remain to screen the wall charge. Such "surface con- 
duction" is driven by the amplified axial electric field in the 
depleted region, which also pushes the negative co-ions to the 
left, thereby sharpening the concentration gradient, leading 
to a steady shock. These effects are illustrated by snapshots 
of (B) counterions and (C) co-ions in a Brownian dynamics 
simulation [36[ . 



II. BASIC PHYSICS OF DESALINATION 
SHOCKS 

Consider the passage of current through a microchan- 
nel with negatively charged side walls, as shown in Fig. [TJ 
Suppose that the EDLs are thin and initially play no role 
in the dynamics. An applied voltage drives current from a 
reservoir on the left to a cation-selective boundary on the 
right, which only allows cations to pass. This boundary, 
shown in Fig. [2j could represent either a cation-selective 
electrodialysis membrane, an electrode where cations are 
reduced to a neutral species, a negative porous electrode 
charging capacitively, or one or more nanochannels with 
over-lapping EDLs. 

In order to maintain electroneutrality as co-ions are ex- 
pelled, the salt concentration is reduced near the bound- 
ary. The ensuing depleted region initially spreads to the 
left by diffusion. As the bulk conductivity is reduced, 
however, the axial electric field is amplified (in order to 
sustain the current) and acts on the counterions screening 
the wall charge to drive surface conduction. Regardless of 
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FIG. 2: Propagation of desalination shock in a straight mi- 
crochannel (A) and a homogeneous microporous medium (B). 
A selective element (membrane) is used at the right-end to 
trigger an initial depletion (concentration polarization effect), 
which then propagates in the form of shock through the mi- 
crostructure. The plot shows the axial profile of the shock 
uniformly sampled in time (C). For a system at constant cur- 
rent, /, and flow rate, Q, the shock propagates at a constant 
speed. 
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the initial Dukhin number, surface conduction eventually 
dominates bulk diffusion in carrying the current through 
the depleted region. Meanwhile, co-ions are driven to 
the left by the large electric field, thus further enhancing 
bulk depletion. This nonlinear feedback causes sharpen- 
ing and propagation of the salt concentration gradient 
similar to standard shock waves. As shown in Fig. [TJA., 
current lines are diverted from the bulk solution into the 
double layers, as they pass through the shock. 

In Fig. [TJ we show the results of Brownian dynamics 
simulations [36|, which clearly illustrate the physics of 
shock propagation. Counterions move from the bulk so- 
lution into the double layers in order to carry current 
around the depleted region behind the shock (Fig. QJ3). 
Meanwhile, co-ions electromigrate ahead of the shock, 
and they become fully depleted behind it (Fig. HP). Al- 
though molecular simulations allow us to visualize the 
trajectories of discrete ions, our goal is to elucidate the 
macroscopic behavior of desalination shocks, so we now 
turn to continuum models. 



tions: 



dci „ „ 
— ^ + u • Vc, = V • 
at 



A(Vc, + pV^) 



(2) 



where we have used the Einstein relation to express the 
mobility of species i as Vi = Di/kT (k = Boltzmann's 
constant, T = absolute temperature) and u is a mean 
fluid velocity in the pores. As a first approximation, 
we have neglected dispersion (velocity-dependent effec- 
tive diffusivity) due to nonuniform convection within the 
pores [45|, Hil , which is reasonable for thin pores [I?} • In 
addition, we enforce macroscopic incompressibility, 

V-u = 0, 

and postulate linear response to gradients of pressure, po- 
tential and concentration at the macroscopic continuum 
scale, 

u = -K H Vp-K E ({ci}, 0) V0 - K D,i({ c i}, <f>) V In c* . 



III. MACROSCOPIC ION TRANSPORT IN 
MICROSTRUCTURES 

The physical arguments above are very general and can 
be extended to microstructures with other geometries. 
As shown in Fig. [2j there is an analogy between macro- 
scopic ion transport in a homogeneous porous medium 
(Fig. [2j>) and in a microchannel (Fig. [2j^) of suitable 
thickness, defined below. We begin by considering uni- 
form microstructures, such as constant-height channels 
and homogeneous porous media (Fig. [2]), and derive gen- 
eral macroscopic transport equations to describe concen- 
tration polarization and desalination shocks. We will 
then extend this model to systems involving geometri- 
cal variations, such as variations in porosity or channel 
cross section. We simply require that the geometrical and 
electrochemical properties of the microstructure vary suf- 
ficiently slowly to justify a volume averaged model. This 
basic assumption also underlies formal homogenization 
analyses [37H42j and leads to macroscopic equations for 
charged porous media of the same general form as we 
propose below [43|, but here we will rely on physical ar- 
guments without deriving any explicit dependence on the 
microstructural geometry. 



The first term is Darcy's law, the second electro-osmotic 
flow, and the third diffusio-osmotic flow, each of which 
in principle have tensorial coefficients in an anisotropic 
medium [48| . The coefficients Ke and Kj^ i depend on 
the ionic concentrations, potential and surface charge and 
could in principle be derived from a microscopic model of 
intrapore transport or approximations for straight chan- 
nels with thin double layers. In our analysis of desalina- 
tion shocks below, we neglect nonlinearities due to con- 
vection to focus on the effects of surface conduction, so 
we leave the derivation and nonlinear analysis of the full 
macroscopic transport equations in three dimensions for 
future work. 



B. Electrostatics 

The key source of nonlinearity in our system is the elec- 
trostatic coupling between ions and the surface charge of 
the microstructure. The electrolyte fills a solid matrix of 
porosity e p (pore volume / total volume) and area den- 
sity a p (pore area / total volume). The walls of the pores 
have a fixed charge density a s (charge / pore area). At 
the macroscopic continuum scale, the surface charge ap- 
pears as a fixed background charge density (charge / pore 
volume) p s given by 



A. Fluxes and flows 

For simplicity, we use dilute solution theory to model 
ionic fluxes, but it is straightforward to extend our re- 
sults by replacing concentrations with activities Q, |44| . 
Let Ci be the mean volume- averaged concentration of ion 
species i of charge qi in the pores (number / pore volume), 
and Di be the effective diffusivity within the porous ma- 
trix [2|, |9|. Conservation of species at the macroscopic 
continuum level is expressed by the Nernst-Planck equa- 



ls = T~ 



(3) 



where h p = e p /a p is an effective pore size. In the first 
step of our derivation, we simply enforce electroneutrality 
at the macroscopic continuum scale, 



e p p + a p a s = 



(4) 



where p is the mean ionic charge density, which is equal 
and opposite to the surface charge density, p s . The 
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macroscopic, volume- averaged electroneutrality condi- 
tion (Eq. 2]) implicitly determines the mean electrostatic 
potential in Eq. [2] This approach has also been employed 
recently to model charge transport in nanochannels [49| 
and carbon nanotubes[50| and can be traced back to early 
models of ion exchange membranes (5l| . 

Let c = \qi\ci be the total ionic charge (regardless 
of sign). For \p s \ <C c, we recover the standard model 
for a quasi- neutral bulk electrolyte, which leads to the 
(ambipolar) diffusion equation for the neutral salt con- 
centration [2|. In the opposite limit, \p s \ ~ c, we recover 
the standard model for a bulk ion-exchange membrane or 
solid electrolytes [521454). In contrast, our focus is on the 
intermediate "leaky membrane" regime, where \p s \ < c, 
which generally introduces nonlinearity due to electromi- 
gration of the diffuse ionic charge that screens the fixed 
background charge. 



C. Binary electrolyte 

We consider the canonical unsupported electrolyte: a 
dilute, asymmetric binary solution (i = +, — ) with arbi- 
trary ionic charges, q± = ±z±e. In this case, macroscopic 
transport equations take the form, 



dc± 



• u • Vc± 




V 2 c±±z±V- (c±V<£) , (5) 



(6) 



where <f> = ecj)/kT is the dimensionless potential, scaled 
to the thermal voltage. Without loss of generality, let us 
assume that the surface charge is negative, p s < 0, and 
use Eq. [6] to replace the ion concentrations c+ and c_ 
with the neutral portion of the salt concentration in the 
bulk (excluding wall shielding charge) 



Cb 



z + c + 



Z-C- 



Ps_ 

e 



2z-c- 



(7) 



In the limit of zero surface charge, this reduces to the 
total concentration of charges (q, — » Z+C+ + z_c_) in a 
neutral electrolyte. In the opposite limit of a fully de- 
pleted bulk electrolyte with nonzero surface charge, this 
quantity vanishes, since only counter- ions remain within 
the EDLs of the microstructure (z + ec + — » —p 3 ). There- 
fore, the variable Cb measures the amount of "free con- 
ductivity" that can be removed from the microstructure 
(i.e. contributing to desalination), without disturbing 
the screening of the fixed surface charge by counter- ions. 
In terms of these variables, the PDEs can be written in 
the following form 



dci) „ — 
-^ + u-Vc 6 = D 
at 







V 2 c 6 - -V • {p a V<fy 
Vj, 



(8) 
(9) 



where j is the volume averaged current density (given 
below); D is the ambipolar diffusivity of a binary elec- 
trolyte [2j (see Appendix [A] for the general form of D and 
z). 



It is clear that in this model, any nonlinear response 
is entirely due to the fixed surface charge, since a linear 
convection- diffusion equation for Cb is recovered from Eq. 
[8] if and only if p s = 0. If any such charge exists in the 
microstructure, then the second term in Eq. [8] survives, 
and the dynamics of the ionic transport will be coupled 
to that of the potential 0, which generally satisfies a PDE 
(Eq. [9]) enforcing the conservation of charge. The nonlin- 
earity becomes apparent from the volume- averaged cur- 
rent density in Eq. [9j which takes the form 



^(j+PsU) 



Kh + T 



V0, 



(10) 



where the second term on the left is the convection of 
charge; the first term on the right is the diffusion current, 
controlled by the parameter 



Da 



z + D + + z-D- 



which measures the asymmetry of the electrolyte; the 
second term on the right hand side of Eq. [10] is Ohm's 
law, where the total conductivity is broken into two parts: 
neutral portion of the bulk, and surface (excess counter- 
ion) contributions. These are respectively: 



(z+i>+ + z-V-)e 2 Cb 



z + v+e\(j s 



(ii) 

(12) 



It is important to stress that what we call n Sl which 
is related to the difference between co- and counter-ion 
concentrations (screening the surface charge), is not the 
same as k' s , the "surface conductivity". The latter is 
defined as the excess conductivity due to sum of co- and 
counter-ion concentrations in the EDLs relative to the 
quasi-neutral bulk solution 0,1110 111. 



IV. CONDUCTIVITY WAVES AT CONSTANT 
CURRENT 

To illustrate the nonlinear dynamics contained in these 
equations, we consider passing a uniform current density 
j = j(t)x and a uniform flow, u = u(t)x through the 
porous medium. We solve Eq. [10] for the electric field and 
substitute back into Eq. [8] to obtain a single, nonlinear 
PDE for bulk conductivity K,b(x,i): 



8Kb 

dt ' 



d_ 

dx 



UKb - 



z-V-e(K, 8 /h p )(j + p s u) 



f^b H~ f^s/hp 



d_ 

dx 



where 

D(K b )=D[l 



D. 



k s j hp 



2z+D+ Kb + Ks/hp 



(13) 



(14) 



This one-dimensional PDE for uniform current is simi- 
lar to that obtained by Mani, Zangle and Santiago [30| 
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in their Simple Model for a flat microchannel with thin 
double layers. Here, we have generalized the model to 
porous microstructures, while adding the convective con- 
tribution of diffuse charge to the current (p s u) as well 
as the conductivity dependence of the effective diffusiv- 
ity D for an asymmetric electrolyte interacting with the 
surface charge. 

If the surface effects (the terms with k s ) can be ne- 
glected, Eq. [13] reduces to the classical linear convection- 
diffusion equation for bulk conductivity. The nonlinear 
flux, Z-V-eK s (j + p s u) I \hpKb-\- k 3 ), can be physically in- 
terpreted as the advection of the surface charge due to 
electromigration (as seen in Eq.[8]). Gradients of this flux 
term are responsible for exchanges between EDL (sur- 
face) and the bulk, which are schematically depicted in 

Fig.m 

Equation [13] has the same form as the equations of 
gas dynamics and shallow water waves [55], and de- 
scribes similar nonlinear wave phenomena. In the long 
time limit in a large system, convection dominates diffu- 
sion and yields a kinematic wave equation of the form, 
Ct + (F(c)) x = 0, which can be solved by the method 
of characteristics. The basic idea is that initial concen- 
tration values propagate with velocity v c = F'(c) along 
characteristic lines in space- time. In order to avoid a 
mult i- valued concentration profile, whenever characteris- 
tics cross, a discontinuity (or shock) in concentration, [c], 
is introduced, which moves at the velocity v s = [F(c)]/[c], 
where [F] is the jump in flux across the shock. The con- 
centration profile across the shock is a traveling wave 
solution, c{x,t) = f(x — v 3 t), to the full equation with 
diffusion. We now apply this kind of analysis to our prob- 
lem. 



A. Dimensionless formulation 

The first step is to define dimensionless variables: 

Kb „ X Z-V-6J ~ t ( Z-V-ej 



D K bo 



D \ K bc 



where fi^oo is the reference bulk conductivity (typically 
in a reservoir connecting to the microstructure). Space 
and time coordinates are nondimensionalized using diffu- 
sive scaling together with characteristic electrodiffusion 
velocity, z-V-ej / ' Kboo- With these definitions, Eq. (p2 
takes the following dimensionless form 



dk 
~dt 



d_ 

dx 



Ps 



d 2 k 
dx 2 ' 



(15) 



where, for simplicity, we have neglected asymmetric dif- 
fusion (D = D) and the convection of diffuse charge 
(\p s u\ <C 1.7 1). In this equation, two fundamental dimen- 
sionless groups appear. The first parameter, 



is the ratio of the mean fluid velocity, u, to the electrodif- 
fusion velocity, Z-V_ej / ' K boo . This parameter affects the 
shock propagation velocity (essentially a Galilean trans- 
formation), but not its dynamics. The second, more im- 
portant, parameter in Eq. [15] is a dimensionless surface 
charge, 



(17) 



hpKjjQ 



h. 



Z- v_ ej ' 



(16) 



With our notation the dimensionless parameter p s in Eq. 
fTTl resembles the Dukhin number, Du, in Eq. [TJ but, as 
discussed above, they are not the same (k s ^ k' s ). 

For typical concentrations in aqueous solutions, p s is 
very small for microstructures (h p ~ 1/im), suggesting 
that the nonlinear term in Eq. [15] can be neglected. One 
mechanism that can activate the nonlinear term (and 
produce shocks) is to locally decrease k to very small 
values of order p s . This is the crucial role that the selec- 
tive boundary (e.g. the membrane in Fig. [2]A.) plays in 
these systems. 

As the shock propagates, it leaves behind a region with 
orders of magnitude lower salt concentration. In other 
words, propagation of the shock acts to desalinate the 
bulk electrolyte. In the next two sections, we analyze 
the dynamics of desalination shocks in systems with non- 
uniform geometries. 



V. WEAKLY VARYING MICROSTRUCTURES 

Figures [3] and [4] show examples of structures involv- 
ing variation of porosity, pore-size, and macroscopic ge- 
ometry. The analysis presented in the previous section 
can be easily extended to these structures. Our analy- 
sis only requires that the microstructure properties vary 
slowly enough to allow a local volume averaged theory. 
While the general derivation is presented in Appendix [A] 
we here continue to focus on the simplified quasi-one- 
dimensional systems and study the response of desalina- 
tion shocks to structural inhomogeneities. Under such 
conditions the modified form of Eq. [13] can be obtained 
by simply scaling all the flux and rate terms with appro- 
priate local volume and/or area measures (see below). 



A. Structures with constant pore size 

We first consider weakly- variable microstructures with 
constant pore size. In other words, in these structures, 
porosity and area-density vary proportionally. With a 
constant surface charge, these structures have a uniform 
background charge density, p s (see Eq. [3]). Figure [3] 
shows examples of such structures, in which the net lo- 
cal volume changes as a function of axial coordinate. In 
Fig. [3]A. the net cross-sectional area (different from area- 
density) is proportional to local porosity, e p ; in Fig.[3£> it 
is proportional to local macroscopic area; and in Fig. [3p 
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Weakly varying microstructures with constant pore 



hp. Schematics include a microstructure with variable 



FIG. 3 

size, 

porosity, e p , and area density, a p , but fixed h p = e p /a p (A), 
a homogeneous microstructure with variation in the macro- 
scopic geometry (B), and a fabricated microchannel with 
variable- width (C). Propagation of a depletion shock through 
the converging- diverging channel under the constant current 
and flow rate condition is shown (D). The plots are sampled 
uniformly in time. 



it is proportional to microchannel width, w. These pa- 
rameters essentially play the same role in modifying the 
dynamics of desalination shocks by scaling the fluxes 
in the conservation laws. For example, for the case of 
the variable-width microchannel, Eq. [13] (again, setting 
p s u = and D = D) will be modified to: 



d , N d 



-V- ejwK s /hp 



dx V dx 



Kb -\- k s I hp 

(18) 

where the "volume averaged" quantities, k^u^ and j are 
effectively the height- averaged quantities , and the equiv- 
alent pore-size, h p , is half the channel height. To be able 
to neglect the transverse fluxes and reduce the system 
to one-dimensional PDE we need the macroscopic geom- 
etry to vary with small slope (dw/dx <C 1), which is a 
standard assumption of lubrication theory. The gradu- 
ally varying assumption imposes an additional condition 
which physically means that macroscopic properties do 
not change much over the axial thickness of the shock. 




FIG. 4: Microstructures with varying pore size, h p . Figure 
shows example of a porous medium (A) and a microtube (B). 
Plots show the bulk conductivity versus axial length across a 
desalination shock (C) for k d = p s = 0.025. Minimum p/a is 
0.25 and is doubled for each subsequent plot up to p/a = 8. 
Figure D shows the plot for p/a = 1 with the dashed lines 
representing the left and right asymptotic curves. 



We use uq : and jo, evaluated at xq (shock location at 
to = 0) to nondimensionalize Eq. [181 w can be nondi- 
mensionalized using wq. Noting that uw and jw are con- 
stant in x due to conservation of mass and current, Eq. 
[18] can be nondimensionalized to 



d_ 

di 



(wk) 



d_ 

dx 



_9_ 

dx 



, dk 
dx 



(19) 



One can verify that Eq. [19] is also applicable to the case 
of porous media. In that case w would be the nondi- 
mensional net cross sectional area. In this formulation u 
and p s are the nondimensional constant parameters, and 
w=w(x) is a known function. Equation fT9l has the trivial 
boundary condition of ft-oo = 1- We also use a Dirichlet 
boundary condition of k{x = 0) = kd = 0(p s ), which 
represents a depletion boundary, initiated by a selective 
element next to the channel. We seek a solution of the 
form 



R(x,t) = f( V ) = f 



X s (t) 



ls(t) 



(20) 



where x s represents the shock location and l s is the shock 
length or axial thickness. The profile of / satisfies an 
ODE, yet to be obtained. Since this profile should look 
like a shock, we have f(rj <C — 1) — 1 ,and f(rj ^> 1) c± kd- 
We propose a solution for x s (t) and l s (t) by speculating 
that the local shock length is proportional to the local 
channel width and its speed is inversely proportional to 
the width: 



dx s v 
di w (x s (t)) 



l s (t) = w (x s (t)) , 



(21) 
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where v is the dimensionless shock speed at x = xq. By 
substituting Eq. [21] into Eq. [20j then into the governing 
equation (Eq. [19]), and ignoring variations of w over the 
shock thickness we obtain the following ODE for /. 



(« - v)f + 



f 



= /"• 



(22) 



To compute the constant v, we can integrate Eq. [22] 
from — oo to +oo and use the boundary conditions. Since 
/' = in the limits, we obtain 



v = u 



0(p s ). 



(23) 



Note that shock propagation would be possible only for 
negative v. This can be typically accommodated only 
if sufficient depletion is introduced at the boundary by 
kd = 0(p s ) (also needs u < 1). 

Substituting Eq. [23] into Eq. [21] and rewriting in the 
dimensional form reveals that for strong shocks (i.e. kd ~ 
p s <C 1) the local shock velocity relative to the local flow 
is 



dx s . v 

" u{x) 



Z-V- ej(x) 



1 



1 + h v Kdjn s 



(24) 



The right-hand-side of Eq. [23] is the electrodiffusion ve- 
locity in the enriched side of the shock scaled by a ra- 
tional function of the surface to bulk conduction in the 
depleted side. As physically expected, in the limit of per- 
fect desalination, kd = 0, the relative shock velocity will 
be identical to the coion electromigration velocity. 
Integrating Eq. [21] yields 



/ 



w(x s )dx s = vt, 



(25) 



which indicates that the rate of sweeping the volume of 
the channel by the shock is constant. This also makes 
sense from the global conservation law point of view: 
Very far from the shock, at the channel boundaries, the 
flux term, uk+{p s )/{k+p 3 ) (see Eq.[T9]h does not change 
with time and the diffusion flux is negligible. From global 
conservation, the depletion of ions inside should balance 
the difference of the fluxes at the boundaries. Therefore, 
the depletion rate should be constant, implying the rate 
of sweeping the volume by the shock should be constant. 



B. Microstructures with variable pore size 

This powerful observation can be generalized to more 
complicated microstructures such as the ones shown in 
Fig. [U In this case, as shown in Fig. HJA., we deal with 
a microstructure with gradually varying porosity, e p and 
surface density, a pi independent of each other. Equiv- 
alently we also can consider microtubal structures (see 
Fig. HP) with gradual variation in cross-sectional area, 



a(x), and cross-sectional perimeter, p(x). Under our sim- 
plifying assumption of quasi-one-dimensional systems, e p 
in the microporous media plays the equivalent role of 
a(x) in microtubal structures; they both scale the bulk 
quantities. In addition, the role of a p in porous media 
is analogous to the role of p(x) in microtubes; they both 
scale the surface quantities. For the case of microtubes, 
the modified governing equation is 



d_ 

dt 



((1Kb) 



d_ 

dx 



Z-V-ejaK 8 /h p \ _ d_ ' ^— dK b 
Kb + Ks/hp ) dx dx 

(26) 

where the "volume averaged" quantities, u, and j are 
effectively the cross-sectional averaged quantities for the 
case of a microtube. The equivalent pore-size, h p , is a/p. 
Equation [26] is very similar to Eq. [18] with the exception 
that now h p is not a constant and is equal to a{x)/p{x). 
Again, as a shock propagates, it sweeps the net available 
volume of the structure at a constant rate independent 
of complexities of a(x) and p(x). 

For the case of constant- we showed that the shock 
axial extent would be proportional to local area of the 
channel. For general a and p however, the evolution of 
shock length is not as simple. It turns out that even a so- 
lution with the form presented by Eq. [20] is not valid any 
more. In this general case, different regions of the shock 
can scale differently. We here only report the analytical 
solution to the shock profile and refer the reader to Ap- 
pendix [B] for details of the derivation. One can show that 
k changes as a function of axial coordinate according to 
the following relation (see Fig. HP) 



x-x s p f M p 

— - — ln(l-«)-(ttd+p s )- In k - n d - 
a a \ a 



Kd 



(27) 

where kd and p s are constants: kd is the dimensionless 
bulk conductivity at the depletion boundary, and p s is 
Ks/hrpix^Kboo- Q> and p are gradually varying local area 
and perimeter nondimensionalized by their reference val- 
ues at Xq. 

With a in the denominator of the left-hand-side, this 
format indicates that the shock axial thickness scales 
with local a (as seen previously), but its shape depends 
on parameter p/a. The right-hand-side of Eq.l27l involves 
two terms: The first term, ln(l — k), is dominant in high 
concentration region (k ^> p s )\ the second term, which 
involves p/a as a parameter, is of order 0(p s ) and is 
dominant in low concentration zone of the shock. A plot 
of the shock profile together with these two asymptotic 
profiles are presented in Fig. [4]C. As mentioned before, 
one can observe that the shock profile is independent of 
convection parameter u. 

From physical standpoint it is worth noting that the 
asymptotic profile of the shock on the high-concentration 
side, 



1 — exp 



Ps 



Kd 



(28) 
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FIG. 5: Schematics of desalination shock propagation in a 
contracting microchannel (7 = —1) is shown in A. Profiles of 
the shock at different stages indicate that as the shock reaches 
the narrower regions of the microchannel it gains speed and 
adopts a sharper axial profile (B). Schematics of propagation 
in a linearly expanding channel is shown in C. Time-series 
of the axial profiles indicate that the shock slows down and 
becomes diffuse towards the end of the channel (D). 



is governed by axial diffusion and a low-concentration 
boundary condition, moving relative to the bulk flow. 
The nonlinear transport associated with surface conduc- 
tivity is negligible through this high-conductivity zone, 
although it plays a role in determining the velocity. The 
same propagating exponential concentration profile of 
Eq. [28] also arises in other situations, such as den- 
dritic electrodeposition (HH, [H3 , where counter- ions are 
removed by convection-diffusion-reaction processes at the 
dendrite tips [58|, rather than by surface conduction. 



VI. SIMILARITY SOLUTIONS FOR 
POWER-LAW GROWTH OF AREA 



TABLE I: Scaling of desalination shock advancement and 
thickening with time for a microchannel with power law 
growth of width. 7 is power of growth of channel width 
with axial coordinate, w = (— x) 7 ; the shock location is as- 
sumed to advance as x s ~ t a ; and the shock axial thickness 
grows/shrinks as l s ~ . 



7 


-1 


(-1,1) 


1 


(1,00) 


a 


exponential 


1 

-y+l 


1 

2 




p 


exponential 


7 
7+1 


1 

2 


1 

2 


description 


shock 


shock 


shock/diffuse 


diffuse 



We seek asymptotic solutions of the form 



(29) 



which describe features that advect with the scaling t a as 
they enlarge (thicken) with the scaling W . Our objective 
is to find a and f3 as functions of 7. Note that a > f3 
would indicate a shock-like solution where propagation 
is faster than growth of the structure; a < f3 indicates 
a diffusion- like spreading, in which advection is not ob- 
servable due to the fast growth of the structure itself. 
Substituting this solution into Eq. [19J and simplifying 
results in 



(Ct a - riP) 1 (Cat a - Prfi p ) + j(Ct a - rfi^)" 



/'+(«/ 
>0) 



In the large t limit appropriately selected a and f3 would 
reduce this equation to an ODE for /. Table U sum- 
marizes the resulting a and /? for different 7 scenarios. 
Following Bazant and Stone [6l|, one can systematically 
check that these are the only scalings that satisfy the 
boundary conditions, but we omit such mathematical 
details here. Note that for the case 7 < —1, the to- 
tal volume of the medium is finite, and an intermediate 
asymptotic limit does not exist. 



A. Intermediate asymptotics 



B. Exponential shock propagation 



In this section we consider the constant-pore-size struc- 
tures again, but with power law growth of their area, 
w — (— x) 7 , as shown in Fig. [U Note that in our nota- 
tion w represents nondimensional cross-sectional area (or 
equivalently channel width or porosity) for a microstruc- 
ture with constant pore size. In this section, variation of 
w is not necessarily negligible over the shock axial extent. 
We are interested in solutions to Eq. [19] at large enough 
times to approach a self-similar form. Such "intermedi- 
ate asymptotic" solutions [59| with power-law monomial 
scalings are expected based on dimensional analysis [60| , 
due to the lack of any natural length scale in the problem. 



In the singular case of 7 = —1 the formal values of 
a and j3 are infinite. Under this condition the correct 
solution would be shock propagation with exponential 
acceleration in time and the correct similarity variable is 
rj = (x + e a t )/e~ a f . In the limit of large t the PDE can 
be transformed to the following ODE: 



(u + a')f + 



f + Ps 



f". 



(31) 



Similar to what observed in Eq. [22j the value of a' can 
be obtained by integrating the above equation from —00 
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to +00 and using the boundary conditions. 
1 



a 



1 + K d /p s 



0(p s ). 



(32) 



The parameter a' can be interpreted as the inverse of the 
time scale for exponential propagation and spreading of 
the concentration profile. 



C. Power- law shock propagation 

For — 1 < 7 < 1 the problem has a power law similarity 
solution with a = 1/(7 + 1) and j3 — 7/(7 + 1). Note 
that for this range a > f3 and thus the solution indicates 
shock propagation. In the limit of large t Eq. [30l reduces 
to the following ODE: 



C 7 + l 

7+1 



f 



(33) 



Interestingly, in the limit of 7 = 1 this solution leads to 
a = /3 = 1/2, which represents the onset of transition 
towards a diffusive propagation. 



D. Diffusive shock propagation in a wedge (critical 

case) 

The case of 7 = 1 represents a structure with linear 
growth of area. A practical example, is a wedge-like chan- 
nel whose width grows with constant slope as shown in 
Fig. [5p. After the case of a straight channel (7 = 0), 
this case maybe the most relevant for lab-on-a-chip sys- 
tems. Note that for 7 = 1 equations can be represented in 
cylindrical coordinates (with x interpreted as radius); the 
lubrication theory assumption (dw/dx <C 1) is not neces- 
sary to enable reduction of the system to one-dimensional 
PDE. Therefore, the wedge angle can be any number 
from to 2tt. 

For 7 = 1 the similarity variable reduces to 77 = xj Vt, 
which shows diffusive scaling in time. Equation [30] re- 
duces to 



This ODE is valid for large t, when the advective flux 
term in Eq. [30] becomes negligible compared to other 
terms. Note that there is no longer any effect of surface 
conduction (p s ) on the intermediate asymptotic similar- 
ity solution. 

In the case that variation of w is due to change in 
the macroscopic geometry, such as in microchannels, for 
very large t the diffusive front may reach locations of 
the channel with large dw/dx and the lubrication theory 
assumption may not be valid any more. As a result Eq. 
l35l will be valid for these structures only for a range in 
time described by: 



1 « « DKb ° 



(36) 



For durations much larger than the upper bound, the 
channel span would have a fast growth, dw/dx ^> 1. In 
this range, the channel maybe approximated by a 180- 
degree wedge and propagation can be modeled by the 
axisymmetric case (7 = 1). 



F. Transients to similarity solutions 

Figure [6] shows a comparison of numerical solutions 
of the full model, Eq. [T9j with our similarity solutions 
for an expanding channel with 7 = 0.5 and a converg- 
ing channel with 7 = —0.25. The spatio temporal plots 
in Fig. [6K and Fig. (6f> show that the shock decelerates 
and becomes smeared by diffusion in the expanding case; 
conversely in the converging channel, the shock sharpens 
and accelerates. Representation of these plots in terms 
of the similarity variable, 77, shows that after a short (di- 
mensionless) transient time the contours collapse into a 
single self-similar profile, as in other problems of inter- 
mediate asymptotics[59|. Comparison with the concen- 
tration profile obtained from the full model demonstrates 
the satisfactory accuracy of the similarity solutions. 



VII. NONLINEAR STABILITY OF 
DESALINATION SHOCKS 



v , i + *A f , _ 1 
2 v J v V/ + p* 



/", (34) 



but there is still some effect of surface conduction, mea- 
sured by p s . 



E. Linear diffusion (no shocks) 

For all values of 7 > 1 the similarity variable will also 
be 77 = x/Vt and Eq. l30l reduces to the following ODE, 
which corresponds to linear diffusion: 



2 + 2),_ r 



(35) 



So far, we have focused on one-dimensional shock pro- 
files, but these are not special cases of the macroscopic 
(volume averaged) nonlinear dynamics. Instead, we ex- 
pect these solutions to be stable attractors, in the sense 
of intermediate asymptotics [59| , at least in the absence 
of flow or sudden property changes (<j s and h p ). To make 
this case, we consider a "thin shock" , whose thickness is 
much smaller than its local radius of curvature, under 
conditions of strong depletion (kd = 0). In this limit, 
the desalinated side contains only surface conductivity, 
thus the Ohm's law in this region would be of the form: 
j = —(K s /h p )V(j). Conservation of charge then implies 
that the potential is harmonic, away from the shock: 



V 2 ^ = for xG Q(t), 



(37) 
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FIG. 6: Spatio-temporal evolution of the desalination shock 
for an expanding microchannel with 7=0.5 (A) and a con- 
tracting microchannel with 7=-0.25 (B). For both channels 
u = 0.5 and p s = 0.1. The black line represents x = 
-ct 1/(7+1) , where c is 0.72 in A and 0.21 in B. When the data 

is plotted against 77 = + c£ 1/(7+1) ) /t 7/(7+1) , the tempo- 
ral evolution collapses to a single profile after sufficient time 
(C,D). Concentration profile at the last time instant (symbol) 
is compared to the asymptotic profile from solution of Eq. [33] 
(E,F). 

where Q(t) represents the desalinated domain. The re- 
gion ahead of the shock has much larger conductivity 
than the desalinated region, so most of the voltage drop 
is sustained in the desalinated region. In this limit, the 
variation of potential outside of Q can be neglected com- 
pared to the scale of potential- variation inside £1: 

= Oforxedfi(*), (38) 

where dfl(t) is the boundary specified by the shock loca- 
tion, x = x s . 

Next, we obtain an equation for boundary- movement 
in terms of potential. As described by Eq.[231 in the limit 
of perfect desalination (k^ = 0), the shock velocity is 
same as local electrodiffusion velocity of the coion species: 

Z-V-e . 

V.s = J- 

^5oo 

Since j is continuous across the shock, it can be written 
using the Ohm's law evaluated at the desalinated side of 
the boundary: 

v. = + fe^W (39) 




FIG. 7: Stability and nonlinear evolution of thin desalination 
shocks in higher dimensions in the absence of flow. The po- 
tential is approximately harmonic in the desalinated region 
behind the shock and constant in the high-conductivity re- 
gion ahead of the shock, and the shock moves in proportion 
to the local electric field, which drives co-ion removal. This 
problem is mathematically equivalent to Laplacian dissolu- 
tion [gj, a well known stable process that leads to smooth 
interfaces from arbitrary initial conditions. 



As shown in Fig. [7J the resulting model is mathemati- 
cally equivalent to the well-known problem of Laplacian 
growth, where an equipotential boundary climbs the nor- 
mal gradient of a harmonic function, only here it is time 
reversed, i.e. the boundary propagates away from the 
harmonic domain. In two dimensions, Laplacian growth 
can be solved using time-dependent conformal maps, and 
it is known to be unstable when the boundary advances 
into the harmonic domain, leading to cusp-like singular- 
ities in finite time [62|. Physically, this situation is like 
dendritic electrodeposition or viscous fingering, where air 
displaces water in a Hele-Shaw cell (without surface ten- 
sion) [63j. In contrast, thin desalination shocks evolve 
by the time-reversed process, which is extremely stable 
and tends to smooth, symmetric shapes. Physically, de- 
salination shock dynamics resemble water displacing air 
in a Hele-Shaw cell or porous medium, or (quasi-steady) 
diffusion-limited dissolution of a porous solid. Dissolu- 
tion fronts are often so stable that they can maintain a 
macroscopic planar shape, even when passing through a 
highly disordered medium [65|, [66|. For several classes 
of analytical solutions of the time-reversed Laplacian 
growth see Ref. [64[. 

This insight justifies a posteriori a key assumption in 
our similarity solutions above. It also shows that they 
represent universally long-time limits for broad classes 
of initial conditions. We leave for future work questions 
of how fluid flow and shock structure might affect this 
picture. Besides microscopic hydrodynamic instabilities 
within the microchannels noted above, we cannot rule out 
the possibility of macroscopic instabilities of desalination 
shocks, e.g. with misaligned fluid flow and electrical cur- 
rent. 
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VIII. CONCLUSION AND OUTLOOK 

In summary, we have developed a general theory of ion 
transport in microchannels and porous media, focusing 
on the new nonlinear regime where surface conduction 
dominates convection in competing with bulk diffusion. 
For slowly varying microstructures, the equations sup- 
port propagating shocks, as well as similarity solutions 
with power-law scalings. Even in the presence of micro- 
scopic inhomogeneities, we expect that these solutions 
are stable attractors of the nonlinear dynamics. The mul- 
tidimensional problem is more complicated, especially in 
situations where the current is misaligned with the fluid 
velocity. We believe this system provides many promis- 
ing directions for research in applied mathematics. 

As suggested by our choice of nomenclature, a natu- 
ral application of our theory would be to water purifi- 
cation and desalination using porous media and mem- 
branes. The basic idea is to extract fresh water continu- 
ously from the region behind a steady desalination shock. 
Our group is currently investigating this concept [67], and 
the results will be reported elsewhere. 

Desalination shocks could also be used to enhance 
the electrokinetic decontamination of microfluidic devices 
and porous rocks, clays or soils [68|,l69|. The propagation 
of a desalination shock would push co-ionic impurities 
ahead of the shock, while counterionic impurities would 
be swept behind the shock by the large electric field. 
This effect, driven by surface conduction, promotes the 
sharpening of the particle profile by electromigration [70| , 
which can also lead to shocks when the particles signifi- 
cantly alter the conductivity [7l|. 

Our theoretical results could also be applied to DC 
electro-osmotic pumps, which employ electro-osmotic 
flow in porous glass frits [72l474j . Strickland et al. [75| and 
Suss et al. [76| have recently found that concentration po- 
larization can be a key factor in the pump performance, 
but current theories do not account for the formation of 
concentration gradients or surface conduction. 

Our results may also find applications in mi- 
cro/nanofluidic systems. We have shown that varying 
the cross-sectional area, perimeter and/or surface charge 
of a microchannel provides robust means to control the 
nonlinear dynamics of transport. In parameter regimes 
where surface conduction is important, this capability 
may be useful in microfluidic devices for biological sam- 
ple pre-concentration [34[ and seawater desalination [35| 
consisting of microchannel/nanochannel junctions. Dur- 
ing normal operation, complex electrokinetic instabilities 
have been observed [15| and, together with fast pressure- 
driven flows [35|, electrohydrodynamic phenomena may 
dominate any effects of surface conduction. Geometrical 
optimization of microchannel interfaces may also lead to 
more robust designs for nanofluidic systems [77j |. e.g. for 
DNA or protein sequencing or molecular sorting, in this 
case to inhibit the formation of shocks, which interfere 
with external control of dynamics within the nanochan- 
nel. 



Another interesting direction would be to relax the as- 
sumption of fixed surface charge, and allow for capaci- 
tive charging [78|, Faradaic reactions [79, 80], or induced- 
charge electro-osmotic flows [8l| in microfluidic devices 
or porous electrodes. Leinweber et al. [82] have observed 
that metal micropost arrays in thin (1 micron) channels 
can produce strong concentration polarization and con- 
tinuous desalination. The effect is driven by surface con- 
duction on ideally polarizable metal cylinders [83]. It is 
likely that desalination shock phenomena, due to surface 
conduction on the microchannel walls, also play a role in 
shaping the salt concentration profile in these devices. 

In the case of porous electrodes, our volume- averaged 
equations for porous media can be applied to capture ef- 
fects of surface conduction, but they must be augmented 
by a charge- voltage relation for the double layer, e.g. us- 
ing the Gouy- Chapman- Stern model of capacitive charg- 
ing |78|,l84| or the Frumkin-Butler-Volmer- Stern model of 
Faradaic reactions [54, 80]. Porous electrodes are widely 
used in electrochemical energy storage devices (batteries, 
super capacitors, fuel cells, etc.) [illzH? but we are n °t 
aware of any prior work considering surface conduction. 
Designing the porous microstructure to exploit the non- 
linear effects of surface conduction could provide a new 
means to enhance the power density of portable power 
sources. 



Appendix A: Porous media with nonuniform 
properties 

In this appendix we present a more general form of 
Eqs.[5J[9j and fTUI applicable to porous media with nonuni- 
form properties such as porosity, diffusivity, and area 
density. We here allow for variable diffusivities, inde- 
pendent of mobility (no Einstein relation). Variable dif- 
fusivity can be due to variable geometrical properties 
of the microstructure or due to nonlinear flow disper- 
sion effects which enhances the effective diffusivity in 
the flow direction [1]. The effect of Taylor dispersion 
due to electro-osmotic flow has been analyzed for thin 
capillaries [85j and flat microchannels [86|, and accu- 
rate volume- averaged equations are available for these 
situations. Yaroschuk and Zholkovskiy [46] have re- 
cently predicted that this effect can also produce sharp 
fronts in the salt concentration in a microchannel, near 
a nanochannel junction, although mainly in thicker mi- 
crochannels (around 100/im) [46]. While the following 
model would accommodate such effects, we here briefly 
note that a simple scaling argument suggests that Taylor 
dispersion can be neglected in very thin (h p <~ /im) mi- 
crostructures due to their relatively low Peclet number, 
Pe = uh p /D 

To derive the model we start with the general form of 
Eq.[5] 



de p c± 
dt 



+V-(e p uc±) = V- 



e p D±Vc± db € p z±v±c±kTV (j) 
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where we remind that e p is the porosity of the porous 
medium. Higher porosity indicates higher effective vol- 
ume to accommodate the transport and thus all fluxes 
scale proportionally with porosity. In this case the con- 
servation laws need to be weighted by the local porosity 
factors. For example, the continuity equation would be 
V • (e p u) = instead of V • u = 0, etc. Rewriting Eq. [A] 
in terms of q>, defined by Eq. [7J and using net neutrality 
(see Eq. [6]) results in 



dc p c b 
dt 



V • (e p uc b ) = V • e p D(Vc h p s V0) + f s 



where, 



= V-(erf), 



— _ Z-V-D+ + z+v+D- 

Z-V- + Z+V+ 



2z+Z-Vj r V-kT 
Z-D+V- + Z+D-V+ 



(Al) 
(A2) 



(A3) 



The f s flux appears as a consequence of nonuniform sur- 
face charge, p s and is equal to 

ep 2£^_ D 

e Z+V+ + Z-V- 

To close the system of Eqs IA1I and IA2I we introduce 
the relation between current and potential gradient, by 
updating Eq. 10 

-j^(}+p s u-D + Vp s ) = -pVn h -[n h + n s /h p ] V0, (A5) 

which only has a slight modification relative to Eq. [10] 
due to nonuniformity of p s with /3, and ft s defined 
the same as in the main text. 



Appendix B: Desalination shock profile in general 
microstructures 



Here we analyze shock structure in a microtubal struc- 
ture whose area a(x) and perimeter p(x) vary indepen- 
dently with position. Due to the mathematical equiva- 
lence of microtubes and porous structures in our model, 
the same analysis also holds for porous medium with vari- 
able porosity e p (x) and surface area density a p (x), which 
respectively play analogous roles as a and p here. We 
start with the nondimensional version of Eq. [26j where 
we use ao and po, respectively the channel cross-sectional 
area and perimeter evaluated at xo , to nondimensionalize 
a and p. 

Using the other dimensionless variables from the main 
text, we arrive at the following dimensionless equation 
describing evolution of bulk conductivity in a channel 
with gradually varying a(x) and p(x): 



where R = Kb/Kboc, % — (%/D)(z-V-ejo/ Kboc), t = 
(t/D)(z-V-ejo/Ki)Oc) 2 , and u = n^n^oc/ Z-V_ejo. To in- 
clude a more general case with gradual variation of sur- 
face conductivity, we define p s = m this 
case p represents variation of both surface charge and 
perimeter and is defined as p = pK s /poK s o. 

We assume that the changes in a and p are slow 
enough, so that their variation over the shock can be 
neglected. We use ki and £2 to denote respectively the 
left and right conductivities out side the shock, but close 
enough so that the cross-section is the same as that at 
the shock. Therefore R\ and R2 may vary as the shock 
sweeps through the channel, which later will be obtained 
from quasi-steady solutions. 

If the shock structure moves with local velocity v, fol- 
lowing the transformation y = x — vt we obtain the fol- 
lowing ODE governing structure of the shock. 



av) + 



VPs 



ah + pp s 



Integration yields 



R{u — av) + — 



VPs 



VPs 



d 

dy 

„ dR 
dy 



, dR 
dy 



C. 



(B2) 



(B3) 



We use Ri and R2 as the boundary condition at infinity. 
Evaluating Eq. IB3I at ±00 and ignoring the diffusion 
term yields the values of C and v: 



(u - aV) 



app s 



(dR 2 +pp s )(dR 1 +pp s )' 



(B4) 



c = vps(dR2 + dRi +yps) , B5 x 

(dR 2 J rpp s )(dR 1 +pp s )' 
Substituting into Eq. IB3I yields: 

dpps {R - k>2)(k - ^1) _ dR 

(dR 2 + pp s )(dRi + pp s ) dR+pp s dy' 

Rearranging terms yields: 

dppgdy dRi + pp s dR dR 2 + pp s dR 

(dR 2 + pp s )(dRi + pps) Ri — R 2 Ri — R R\ — R 2 R — R 2 

(B7) 

Integration results in 

dpp s (y-y ) dRi+pps dR 2 + pp s 

-jz-i w _ t^t — ~ — ln(«i-«) ~ ~ — ln(«- 

[an 2 +pp s )(aKi + pp s ) k\ — k 2 k\-k 2 

(B8) 

Now we need to substitute values of R\ and R 2 in terms Rd 
and local p and a. R 2 satisfies the steady state condition 
for Eq. IB1I in the depletion region. Since we are far 
from the shock the diffusive flux can be neglected in this 
region; hence the net convective flux should be constant 
in order to satisfy the steady state condition. Therefore, 
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Note that p and a are one at x = xq. Considering the 
fact that k>2 ~ K>d ~ 0{ps) <C 1, we can simplify this 
expression and arrive at 



= + o{p z s ). 

Similarly, one can show that 

«i = l + 0(p a ). 



(BIO) 



(Bll) 



Substituting these expressions for k\ and into Eq. 
results in 



y-Vs 



ln(l-k)-(k d +p s )(p/a)ln [k - k d {p/a] 

K Kd + Ps J a 

(B12) 

which is a direct relation between the bulk conductivity 
and axial coordinate across a shock. Having x s = yo + vt 
this equation can be transformed to Eq. [271 

Figure HP shows the shock profiles obtained from Eq. 
IB12I One can see that different regions of the shock scale 
differently as parameters a and p vary. While the high- 
concentration region of the shock scales with local a, the 
low-concentration region is dependent on both parame- 
ters a and p. This also makes sense from the form of Eq. 



IB12l since the high- and low-concentration regions can be 
approximated respectively by the first and second term 
in the right hand side of the Eq. IB12I A plot of the 
shock profile together with these two asymptotic profiles 
are shown in Fig. |4j 

In practical scenarios the conductivity-drop across the 
shock is orders of magnitude (0(p s ) <C 1). Under such 
conditions most of the drop, from £ = ltop s <C£<Cl, 
can be approximated by only the first term on the right- 
hand-side of Eq. IB12[ Therefore, as a rule of thumb, one 
can say that the shock thickness approximately scales 
with local area. Note that this simple criterion assumes 
that variations in p/a are finite and bounded with an 
] upper bound much smaller than l/p s - 
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